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Abstract 

We develop a phase-field model of eutectic growth that uses three phase 
fields, admits strictly binary interfaces as stable solutions, and has a smooth 
free energy functional. We use this model to simulate oscillatory limit cycles 
in two-dimensional lamellar growth, and find a continuous evolution from low- 
amplitude oscillations to successive invasions of one solid phase by the other 
when the lamellar spacing is varied. 

I. INTRODUCTION 

Solidification is both a fascinating example of pattern formation out of equilibrium and 
a phenomenon of practical importance in metallurgy. The most common solidification mi- 
crostructures found in industrial alloys are dendrites and eutectic composite structures (rods 
or lamellae). Eutectic growth occurs when two solid phases of different composition can solid- 
ify from the same melt. The interplay between the redistribution of chemical components in 
front of the moving phase boundary and the effects of capillarity along the curved interfaces 
gives rise to a wealth of different patterns and nonlinear phenomena such as bifurcations, 
limit cycles, solitary waves, and even chaotic states JT|]. 

Phase-field modeling has become the method of choice for simulating solidification fronts. 
Its main advantage is that it avoids explicit tracking of the solidification front by introduc- 
ing "phase fields" that vary continuously between constant values corresponding to the bulk 
phases, thus replacing the sharp fronts by diffuse interfaces with a finite thickness. The 
connection to the traditional free-boundary formulation of solidification is established by 
the technique of matched asymptotic expansions around the equilibrium front in the limit 
where the thickness of the diffuse interfaces is much smaller than all other relevant length 
scales. Recently, the computational efficiency of phase-field models for single-phase solidifi- 
cation of pure substances and binary alloys has been drastically enhanced by pushing this 
perturbation analysis to second order ||. This, combined with a random walk algorithm 
for an efficient simulation of the diffusion equation, has made it possible to simulate quan- 
titatively the dendritic growth of a pure substance in three dimensions for experimentally 
relevant parameters ||. 

Therefore, it would be highly desirable to extend this second-order asymptotic analysis 
to phase-field models for multi-phase solidification; we will specifically address the case of 
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eutectic growth, where two solid phases denoted by a and (3 grow from the liquid. The first 
phase-field models for eutectic growth used the standard phase field to distinguish between 
solid and liquid, and a concentration field j|,|5| or a second phase field || to distinguish 
between the two solids. However, a solid-liquid interface then involves variations of the two 
variables and the asymptotic analysis becomes extremely cumbersome. The more recent 
multi-phase-field approach |7j] assigns one phase field to each phase and interprets these 
fields as local volume fractions. Across an interface, in principle only the volume fractions of 
the two bulk phases limiting the interface need to vary, and one of them can be eliminated in 
terms of the other. Hence, this method allows for equilibrium interface solutions depending 
on a single variable that are equivalent to the usual "binary" interfaces of the standard 
phase-field model. However, to achieve this the model has to be carefully designed, and so 
far strictly binary interfaces have only been obtained by using a singular free energy (double 
obstacle potential) ||, which might complicate the asymptotic analysis. 

We present a phase-field model for eutectic growth that uses three phase fields, yields 
exactly binary interfaces at equilibrium, and has a smooth free energy functional. Fur- 
thermore, the free energy is designed to keep the interfaces binary to first order in the 
usual non-equilibrium asymptotic expansion, which makes it a promising starting point for 
a second-order analysis. 

To illustrate the capabilities of the model in its present form, we simulate lamellar eutectic 
growth of a simple model alloy with a symmetric phase diagram. We show that small- 
amplitude oscillations of the lamellar width, "giant oscillations" , and successive invasions of 
one solid phase by the other, each of which is seen in experiments, all belong to the same 
branch of oscillatory limit cycles that is parametrized by the lamellar spacing. 



II. MODEL 



We use three phase fields pi G [0, 1] representing the local volume fractions of each phase 
(a, (3 or liquid), and thus J2iPi — 1- Their dynamics is derived from a free energy functional 
J 7 , 



dpi 
dt 
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(1) 



where r is a relaxation time and the variational derivative has to take into account the 
constraint Y^,iPi — 1- Furthermore, the concentration field c obeys the conserved dynamics 



dc/dt = V M(p a , PfB , Pl )Vii 



(2) 



where /i = 5 J- '/5c is the chemical potential and M is a phase-dependent mobility. 

The free energy functional J- = J v fdV is the volume integral of a free energy density 
/ that is conveniently decomposed into / = / gra d + /tw + fc, where /tw is a triple well 
potential with minima on each pure phase — the equivalent of the double well potential of 
the standard phase-field model, / gra d contains the gradients of the phase fields and ensures 
a finite interface thickness and surface tension, and f c couples the phase fields and the 
concentration and drives the system out of equilibrium. 
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FIG. 1. Triple- well potential fxw ( a ) an d elementary tilt function —gi (b) drawn as "land- 
scapes" over the Gibbs simplex. 

Due to the constraint J2iPi — 1; the f ree energy density and the interface solutions can 
be conveniently represented on a Gibbs simplex, that is an equilateral triangle where the 
distance to each side of the triangle from a given point represents the value of one of the 
phase fields. Thus, each vertex corresponds to a pure phase, and each side to a purely binary 
interface. Our goal is to obtain strictly binary interfaces, i.e., the equilibrium solutions of 
(|I|) for an interface connecting two phases should exactly run along the edges of the simplex. 
One way to achieve this is to choose 



and to design the function fxw m such a way that the three minima are connected by 
"saddles" that run along the edges and that have vanishing derivatives in the direction 
normal to the edges. Such functions can be constructed using the geometry of the Gibbs 
simplex, as will be detailed elsewhere. The simplest choice, 



is plotted in Fig. |T]a and reduces to the standard fourth-degree double well potential of the 
binary phase-field model on the edges. With these choices, the free energy is completely 
symmetric with respect to the interchange of any two phases, and as a consequence all three 
surface tensions are equal. However, the general case of unequal surface tensions can also 
be treated by adding other terms in the free energy density that modify the height of the 
saddle points. 

Out of equilibrium, the triple well is "tilted" by the term f c in order to favor one phase 
over the others. Since we want to keep the interfaces binary, this term needs to be carefully 
constructed. The elementary building blocks are tilt functions that satisfy gi — 1 for 
Pi = 1, gi = for pj = 1, j ^ i, and that have vanishing derivatives in the direction normal 
to the edges of the simplex. One possible choice is 



grad 




(3) 




(4) 
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;i + Pj - Pi) 2 (2 + Pi - pj) 2 + (i + Pk - Pi) 2 (2 + Pl ~ Pk f 

+ (Pk~P,) 2 (3-\p k -p,\) 2 ] /12 . 



(5) 



In Fig. |l]b, we plot —gi that lowers the well for the liquid with respect to the two others. 
With the help of these functions, 



fc 



(6) 



interpolates between the concentration-dependent free energies fi = (c — Ai) 2 /2 + Bi of each 
phase i. Here, c = (C — C#)/AC is a scaled concentration variable, where AC = Cp — C a 
is the width of the eutectic plateau, and C a , Cp, and Ce are the equilibrium compositions 
of the two solids and the liquid at the eutectic temperature. An arbitrary phase diagram 
can be constructed by performing the well-known double tangent construction, which yields 
cf = Ai + (Bj — Bi)/(Aj — A^ for the equilibrium concentration cf of phase i coexisting 
with phase j, and by choosing appropriate Aj(T), B^T), where T is the temperature. Note 
that by the same procedure, a peritectic phase diagram can be constructed ||. 

Finally, we need to specify the mobility function in (0). Since d 2 fi/dc 2 = 1, the choice 
M(p a ,pp,pi) = Dpi yields a constant diffusivity D in the liquid, and no diffusion in the solid 
(one-sided model). 



III. SIMULATIONS 

In thin-sample directional solidification of eutectic alloys with global composition in 
the eutectic range, the basic pattern is a periodic array of alternating lamellae of each 
solid phase, growing perpendicular to the large-scale solidification front and parallel to the 
imposed temperature gradient. This basic state exhibits various instabilities that have been 
carefully studied both experimentally in the transparent organic alloy CBr 4 -C 2 Cl 6 pj and 



by numerical simulations using the boundary integral method |]10 |. For fixed experimental 
parameters, the stability of the basic state is controlled by the width of one lamella pair, or 
lamellar spacing. For sample compositions close to the eutectic point, the first instability 
that is encountered for increasing spacing is a period-preserving oscillatory instability, that 
is, the lamellae start to oscillate, with all the lamellae of the same phase oscillating in phase, 
and the global spacing left unchanged. Beyond the onset of this instability, stable limit 
cycles are found, with an oscillation amplitude that increases with spacing. In the boundary 
integral simulations, this branch of solutions terminates when the amplitude becomes too 
large because the thinner lamellae pinch off. In contrast, in experiments "giant oscillations" 
with very large amplitudes can be obtained [|ll|]. Furthermore, when the experiments are 



started from a single solid phase growing into the liquid with a few widely spaced nuclei of 
the other solid phase on the solidification front, the nuclei grow and spread along the front 
forming "invasion tongues" ||12|| . When two tongues growing in opposite directions collide, 
they do not coalesce; instead, a narrow channel of the other phase remains and the whole 
process starts over with the role of the two phases reversed. The interface dynamics is hence 
still oscillatory, even if the patterns look completely different. In this regime, the boundary 
integral method becomes inapplicable because it uses the quasistationary approximation of 
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the diffusion equation that is valid only if the interface motion is slow compared to the 
diffusive solute redistribution. Here, we investigate oscillatory limit cycles with our phase- 
field model, that does not have this restriction. 

The sample is solidified with pulling speed V in a constant temperature gradient G 
directed along the z direction. For simplicity, we use a symmetric phase diagram that 
has constant concentration jumps across the interfaces (parallel liquidus and solidus lines). 
This is implemented by choosing A a = —0.5, Ap = 0.5, A t = 0, B a = Bp = and 
Bi = — (z — Vt)/l T , where It = (mAC)/G is the thermal length, with m being the liquidus 
slope (equal for both phases). The other relevant physical length scales are the diffusion 
length = D/V and the capillary length do = r/(mAC), where T is the Gibbs-Thomson 
constant. In terms of the model parameters, d = (2\/2/3)W. All simulations discussed 
here were performed with r = W — D — 1, It /do = 530, Id /do = 212 and at exactly eutectic 
composition (that is, c = far ahead of the front). 

The equations are integrated by a standard explicit finite-difference scheme on a regular 
grid of spacing 0.8 W for which we have checked that the discretization has converged. 
Since the diffusivity vanishes in the solid, the diffusion equation only needs to be solved in 
interfacial and liquid regions. Moreover, far ahead of the solidification front, it is solved 
using a variant of the random walker algorithm of Ref. || that will be detailed elsewhere. 
Two types of initial conditions are used. To obtain lamellar steady states, simulations are 
started from two flat lamellae in contact with the liquid. If the lamellar state is unstable, 
the numerical noise generated by the random walkers triggers the instability. Alternatively, 
simulations are started from a flat single solid phase at equilibrium with the liquid. Upon 
pulling, the interface recoils towards colder temperatures as a diffusion layer of the rejected 
component is gradually formed. When the undercooling reaches a predetermined value, 
nucleation is mimicked by placing a large enough semicircular nucleus of the other phase on 
top of the interface. 

We restrict our attention here to period-preserving oscillations without lateral drift. 
Therefore, it is sufficient to simulate two adjacent half lamellae or half a nucleus of one 
solid phase on top of the other with reflecting (no-flux) boundary conditions on the sides 
of the simulation box that are parallel to the lamellae. The lamellar spacing A is varied by 
changing the lateral size of the simulation box. In the figures below, for clarity we have 
reconstructed a whole lamella pair by reflecting the system with respect to one side of the 
box. 

We start by constructing the branch of stable lamellar steady-state solutions and find 
the minimal undercooling spacing at A m i n ~ 85<io, in good agreement with the prediction 
of the Jackson-Hunt theory [|T^], A m i n ~ 79d . Next, we start from a single nucleus and 
monitor the dynamics for increasing spacings. We find a bifurcation toward oscillatory limit 
cycles for A close to 2A m i n . As shown in Fig. ||a, for A = 2A m i n the oscillation already has 
an amplitude of about a quarter of the lamellar spacing, but its shape is still close to a 
sinusoidal wave. At A = 4A m ; n (Fig. 2 [|b), the oscillation is highly asymmetric, and has an 
amplitude close to half the lamellar spacing, i.e., only a thin channel of the other solid phase 
is left. Finally, at 8.48A min (Fig. ||c), three different regimes become apparent: (i) initially, 
the nucleus grows slowly and spreads along the interface to form an invading finger with a 
well-defined shape, (ii) the finger speeds up with an approximately constant acceleration, 
leaving a straight border with the other solid phase behind, and (iii) when it approaches a 
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FIG. 2. Snapshot pictures of typical oscillatory structures. Growth is upwards, the two solids 
are black and white, and the greyscale in the liquid is proportional to the solute concentration, 
(a) low-amplitude oscillations (A/A m i n = 2), (b) giant oscillations (A/A m i n = 4), (c) successive 
invasions (A/A m i n = 8.48). The scale is the same on all three figures. 



finger growing in the opposite direction (in our simulations, its mirror image generated by 
the reflecting boundary conditions), it "feels" the diffusion field generated by the other and 
slows down very rapidly, allowing for the other phase to emerge through a narrow channel 
left between the two fingers and thus restart the process with the role of the two solids 
reversed. 

Remarkably, our largest simulations reproduce many characteristic features of the "suc- 
cessive invasions" observed in the experiments ||12|| , in particular the approximately constant 
acceleration of invasion tongues and the leftover of narrow channels between "colliding" fin- 
gers. From our simulations, it hence appears that this regime can be understood as period- 
preserving oscillations with very large amplitudes. Indeed, the morphology at 4A min already 
presents stages (i) and (iii), but not (ii). Apparently, to observe this intermediate, constant 
acceleration regime, it is sufficient to have enough "free space" left in front of the propa- 
gating finger (i.e., large lamellar spacings). It is also interesting to note that in peritectic 
alloys, simulations of spreading in the same geometry produce stages (i) and (ii), but then 
colliding fingers coalesce, and to obtain an oscillatory dynamics explicit nucleation events 
have to be introduced PJI3|. 

The fact that we can switch from gentle oscillations to long invading fingers just by 
changing the lamellar spacing demonstrates that both belong to the same branch of oscilla- 
tory limit cycles. This branch of solutions is indeed an attractor, since simulations at fixed 
spacing and started from very different initial conditions always converge to the same cycle. 



IV. CONCLUSION 



We have developed a phase-field model of eutectic solidification that has a smooth free 
energy landscape and yields exactly binary interfaces away from the trijunction points. 
These properties make it a promising starting point for a second-order asymptotic analysis. 
Furthermore, we have used the model to simulate oscillatory limit cycles in a model alloy 
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that has a symmetric phase diagram, and we have found that low-amplitude oscillations, 
"giant oscillations", and successive invasions, all observed in experiments, lie on a single 
branch of solutions that is parametrized by the lamellar spacing and that bifurcates from 
the steady-state branch at about twice the minimum undercooling spacing. 

Here, we have focused on period-preserving oscillatory modes; many other instability 
modes exist, for example the tilt instability or period-doubling oscillations. These insta- 
bilities can be studied using less restrictive boundary conditions. Preliminary simulations 
show that for the case studied here, namely an alloy with symmetric phase diagram at its 
eutectic composition, period-preserving oscillations remain the only stable limit cycles. In 
contrast, if we use off-eutectic sample compositions or an asymmetric phase diagram, and 
hence break the complete symmetry between the two solid phases, other modes become 
active. Therefore, our model can be used to obtain a more complete bifurcation diagram, 
which is the subject of ongoing work. 

We thank Silvere Akamatsu and Gabriel Faivre for many stimulating discussions and 
Jean-Frangois Gouyet for his help with Fig. |I|. R. F. acknowledges financial support from 
project FMRX-CT96-0085 of the European Commission and from Centre National d'Etudes 
Spatiales (France). 
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